# Supplementary R script
# Multiple Correspondence Analysis (MCA) of climate-related stressors,
# ecosystems, locations and risk types
#
# Purpose:
#   This script reproduces two MCA analyses used to explore associations among
#   categorical variables in the reviewed literature.
#
# Input files:
#   1) MCA.xlsx      - first categorical matrix
#   2) ACM 2.xlsx    - second categorical matrix
#
# Output files:
#   1) MCA_results_Prueba1_standardised_no_rare_no_reported_corrected.xlsx
#   2) MCA_categories_Prueba1_standardised_no_rare_no_reported_Benzecri_corrected.tiff
#   3) MCA_results_Prueba2_standardised_no_rare_no_reported_corrected.xlsx
#   4) MCA_categories_Prueba2_standardised_no_rare_no_reported_Benzecri_corrected.tiff
#
# Notes:
#   - Categories labelled as "Not reported" or equivalent terms are removed.
#   - Rare categories are removed when they occur fewer than min_n times.
#   - MCA axes in the figures are labelled using Benzécri-corrected inertia.
#   - Greenacre values are reported here as percentages relative to the total
#     positive corrected inertia, following the same implementation used for
#     the Benzécri-corrected percentages in this script.
#
# Required R packages:
#   readxl, dplyr, stringr, FactoMineR, openxlsx, ggplot2, ggrepel, tidyr
# ============================================================


# ============================================================
# 1. USER SETTINGS
# ============================================================

# Set the folder containing the input Excel files.
# Change this path before running the script on another computer.
work_dir <- "C:/Users/xx/xx"
# Input file names and sheet names
file_matrix_1 <- "MCA.xlsx"
file_matrix_2 <- "ACM 2.xlsx"
sheet_name <- "Hoja1"

# Category filtering
remove_rare_categories <- TRUE
min_n <- 2

# Figure export settings
figure_dpi <- 300
figure_compression <- "lzw"


# ============================================================
# 2. PACKAGES
# ============================================================

packages <- c(
  "readxl", "dplyr", "stringr", "FactoMineR",
  "openxlsx", "ggplot2", "ggrepel", "tidyr"
)

packages_to_install <- packages[!packages %in% installed.packages()[, "Package"]]

if(length(packages_to_install) > 0){
  install.packages(packages_to_install)
}

library(readxl)
library(dplyr)
library(stringr)
library(FactoMineR)
library(openxlsx)
library(ggplot2)
library(ggrepel)
library(tidyr)

setwd(work_dir)


# ============================================================
# 3. GENERAL FUNCTIONS
# ============================================================

# Clean text strings by removing non-breaking spaces, replacing different
# dash types with a standard hyphen, and trimming repeated spaces.
clean_text <- function(x){
  x %>%
    as.character() %>%
    str_replace_all("\u00A0", " ") %>%
    str_replace_all("\u2011", "-") %>%
    str_replace_all("\u2013", "-") %>%
    str_replace_all("\u2014", "-") %>%
    str_squish()
}


# Standardise column names by removing non-breaking spaces and repeated spaces.
clean_column_names <- function(x){
  x %>%
    str_replace_all("\u00A0", " ") %>%
    str_squish()
}


# Replace categories occurring fewer than min_n times with NA.
remove_rare_levels <- function(x, min_n = 2){
  tab <- table(x)
  rare_levels <- names(tab[tab < min_n])
  ifelse(x %in% rare_levels, NA, x)
}


# Calculate raw inertia and corrected inertia following the implementation
# used in this analysis.
mca_corrected_inertia <- function(mca_res, data_mca){

  Q <- ncol(data_mca)

  eig_raw <- as.data.frame(mca_res$eig)
  eig_raw$Dimension <- rownames(eig_raw)
  rownames(eig_raw) <- NULL

  lambda <- eig_raw$eigenvalue
  threshold <- 1 / Q

  # Benzécri correction
  benzecri_eigen <- rep(0, length(lambda))
  valid <- lambda > threshold

  benzecri_eigen[valid] <- (Q / (Q - 1))^2 *
    (lambda[valid] - threshold)^2

  benzecri_percent <- rep(0, length(lambda))

  if(sum(benzecri_eigen) > 0){
    benzecri_percent <- benzecri_eigen / sum(benzecri_eigen) * 100
  }

  # Greenacre values are reported as percentages relative to the total
  # positive corrected inertia, matching the analytical implementation
  # used in this script.
  greenacre_eigen <- benzecri_eigen

  greenacre_percent <- rep(0, length(lambda))

  if(sum(greenacre_eigen) > 0){
    greenacre_percent <- greenacre_eigen / sum(greenacre_eigen) * 100
  }

  data.frame(
    Dimension = eig_raw$Dimension,

    Raw_eigenvalue = lambda,
    Raw_percent = eig_raw$`percentage of variance`,
    Raw_cumulative = eig_raw$`cumulative percentage of variance`,

    Benzecri_eigenvalue = benzecri_eigen,
    Benzecri_percent = benzecri_percent,
    Benzecri_cumulative = cumsum(benzecri_percent),

    Greenacre_eigenvalue = greenacre_eigen,
    Greenacre_percent = greenacre_percent,
    Greenacre_cumulative = cumsum(greenacre_percent)
  )
}


# Export the main MCA outputs to an Excel workbook.
export_mca_results <- function(
    output_file,
    data_mca,
    mca_res,
    corrected_inertia
){

  eig <- as.data.frame(mca_res$eig)
  eig$Dimension <- rownames(eig)
  rownames(eig) <- NULL
  eig <- eig[, c("Dimension", colnames(eig)[1:3])]

  category_contributions <- as.data.frame(mca_res$var$contrib)
  category_contributions$Category <- rownames(category_contributions)
  rownames(category_contributions) <- NULL

  category_cos2 <- as.data.frame(mca_res$var$cos2)
  category_cos2$Category <- rownames(category_cos2)
  rownames(category_cos2) <- NULL

  category_coordinates <- as.data.frame(mca_res$var$coord)
  category_coordinates$Category <- rownames(category_coordinates)
  rownames(category_coordinates) <- NULL

  individual_coordinates <- as.data.frame(mca_res$ind$coord)
  individual_coordinates$Article_ID <- 1:nrow(individual_coordinates)

  wb <- createWorkbook()

  addWorksheet(wb, "Clean_standardised_data")
  writeData(wb, "Clean_standardised_data", data_mca)

  addWorksheet(wb, "Eigenvalues_raw")
  writeData(wb, "Eigenvalues_raw", eig)

  addWorksheet(wb, "Corrected_inertia")
  writeData(wb, "Corrected_inertia", corrected_inertia)

  addWorksheet(wb, "Contributions")
  writeData(wb, "Contributions", category_contributions)

  addWorksheet(wb, "Cos2")
  writeData(wb, "Cos2", category_cos2)

  addWorksheet(wb, "Coordinates_categories")
  writeData(wb, "Coordinates_categories", category_coordinates)

  addWorksheet(wb, "Coordinates_articles")
  writeData(wb, "Coordinates_articles", individual_coordinates)

  saveWorkbook(wb, output_file, overwrite = TRUE)
}


# ============================================================
# 4. MCA ANALYSIS - MATRIX 1
# ============================================================

# Matrix 1 variables:
#   Primary_stressor: main climate-related stressor category
#   Ecosystem_type: ecosystem type
#   Co_stressor: non-climatic or additional stressor category
#   Risk_Type: type of risk

data_1 <- read_excel(file_matrix_1, sheet = sheet_name)

names(data_1) <- clean_column_names(names(data_1))

data_mca_1 <- data_1 %>%
  select(
    Primary_stressor,
    Ecosystem_type,
    Co_stressor,
    Risk_Type
  ) %>%
  mutate(across(everything(), clean_text)) %>%
  filter(if_all(everything(), ~ !is.na(.) & . != ""))


# ------------------------------------------------------------
# 4.1 Standardise categories - Matrix 1
# ------------------------------------------------------------

data_mca_1 <- data_mca_1 %>%
  mutate(

    Primary_stressor = case_when(
      Primary_stressor %in% c(
        "Extreme weather event",
        "Extreme weather events",
        "Extreme events"
      ) ~ "Extreme weather events",

      Primary_stressor %in% c(
        "Marine Heatwaves",
        "Marine heatwaves",
        "Marine heatwave"
      ) ~ "Marine heatwaves",

      Primary_stressor %in% c(
        "Sea level rise",
        "Sea-level rise",
        "SLR"
      ) ~ "Sea level rise",

      Primary_stressor %in% c(
        "Climate change",
        "Climate variability",
        "Climate change and variability",
        "Climate variability and change",
        "Climate change / variability"
      ) ~ "Climate change / variability",

      TRUE ~ Primary_stressor
    ),

    Ecosystem_type = case_when(
      Ecosystem_type %in% c("Coral reef", "Coral reefs") ~ "Coral reefs",
      Ecosystem_type %in% c("Mangrove", "Mangroves") ~ "Mangroves",
      Ecosystem_type %in% c("Kelp forest", "Kelp forests") ~ "Kelp forests",
      Ecosystem_type %in% c(
        "Shelf and pelagic",
        "Shelf/pelagic",
        "Shelf & pelagic"
      ) ~ "Shelf and pelagic",
      Ecosystem_type %in% c("Sandy beach", "Sandy beaches") ~ "Sandy beaches",
      TRUE ~ Ecosystem_type
    ),

    Co_stressor = case_when(
      Co_stressor %in% c(
        "No reported",
        "Not reported",
        "None reported",
        "Not specified",
        "NA"
      ) ~ "Not reported",

      Co_stressor %in% c(
        "Natural / physical processes",
        "Natural/physical processes",
        "Natural physical processes"
      ) ~ "Natural/physical processes",

      Co_stressor %in% c(
        "Land use",
        "Land-use",
        "Land use change",
        "Land-use change"
      ) ~ "Land use",

      Co_stressor %in% c(
        "Urbanisation",
        "Urbanization"
      ) ~ "Urbanisation",

      TRUE ~ Co_stressor
    ),

    Risk_Type = case_when(
      Risk_Type %in% c(
        "Socio ecological",
        "Socio-ecological",
        "Socioecological"
      ) ~ "Socio-ecological",

      Risk_Type %in% c(
        "Biogeophysical",
        "Bio-geophysical"
      ) ~ "Biogeophysical",

      TRUE ~ Risk_Type
    )
  )


# ------------------------------------------------------------
# 4.2 Remove not reported values and rare categories - Matrix 1
# ------------------------------------------------------------

data_mca_1 <- data_mca_1 %>%
  filter(
    Primary_stressor != "Not reported",
    Ecosystem_type != "Not reported",
    Co_stressor != "Not reported",
    Risk_Type != "Not reported"
  )

if(remove_rare_categories){
  data_mca_1 <- data_mca_1 %>%
    mutate(across(
      everything(),
      ~ remove_rare_levels(., min_n = min_n)
    )) %>%
    drop_na()
}

data_mca_1 <- data_mca_1 %>%
  mutate(across(everything(), as.factor))

print(lapply(data_mca_1, table))
str(data_mca_1)
summary(data_mca_1)


# ------------------------------------------------------------
# 4.3 Run MCA - Matrix 1
# ------------------------------------------------------------

mca_1 <- MCA(data_mca_1, graph = FALSE)

corrected_inertia_1 <- mca_corrected_inertia(mca_1, data_mca_1)

print(corrected_inertia_1)


# ------------------------------------------------------------
# 4.4 Export results - Matrix 1
# ------------------------------------------------------------

export_mca_results(
  output_file = "MCA_results_Prueba1_standardised_no_rare_no_reported_corrected.xlsx",
  data_mca = data_mca_1,
  mca_res = mca_1,
  corrected_inertia = corrected_inertia_1
)


# ------------------------------------------------------------
# 4.5 Prepare plot data - Matrix 1
# ------------------------------------------------------------

coord_plot_1 <- as.data.frame(mca_1$var$coord)
names(coord_plot_1)[1:2] <- c("Dim1", "Dim2")
coord_plot_1$Category_full <- rownames(coord_plot_1)

# FactoMineR may return category names either with or without variable prefixes.
# The following mapping handles both possibilities.
map_prefixed_1 <- bind_rows(
  lapply(names(data_mca_1), function(v){
    data.frame(
      Category_full = paste0(v, "_", levels(data_mca_1[[v]])),
      Category = levels(data_mca_1[[v]]),
      Variable_original = v
    )
  })
)

map_raw_1 <- bind_rows(
  lapply(names(data_mca_1), function(v){
    data.frame(
      Category_full = levels(data_mca_1[[v]]),
      Category = levels(data_mca_1[[v]]),
      Variable_original = v
    )
  })
) %>%
  group_by(Category_full) %>%
  filter(n() == 1) %>%
  ungroup()

map_variables_1 <- bind_rows(map_prefixed_1, map_raw_1)

coord_plot_1 <- coord_plot_1 %>%
  left_join(map_variables_1, by = "Category_full") %>%
  mutate(
    Variable = recode(
      Variable_original,
      Primary_stressor = "Primary stressor",
      Ecosystem_type = "Ecosystem type",
      Co_stressor = "Synergistic stressors",
      Risk_Type = "Risk type"
    )
  )

coord_plot_1$Variable <- factor(
  coord_plot_1$Variable,
  levels = c(
    "Primary stressor",
    "Ecosystem type",
    "Synergistic stressors",
    "Risk type"
  )
)

table(coord_plot_1$Variable, useNA = "ifany")

dim1_1 <- round(corrected_inertia_1$Benzecri_percent[1], 1)
dim2_1 <- round(corrected_inertia_1$Benzecri_percent[2], 1)


# ------------------------------------------------------------
# 4.6 Plot and export figure - Matrix 1
# ------------------------------------------------------------

p_mca_1 <- ggplot(
  coord_plot_1,
  aes(x = Dim1, y = Dim2, colour = Variable)
) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.4) +
  geom_point(size = 1.8, alpha = 0.8) +
  geom_text_repel(
    aes(label = Category),
    size = 3,
    max.overlaps = Inf,
    box.padding = 0.25,
    point.padding = 0.15,
    segment.size = 0.2,
    show.legend = FALSE
  ) +
  scale_colour_manual(values = c(
    "Primary stressor" = "#6A5ACD",
    "Ecosystem type" = "#4D4D4D",
    "Synergistic stressors" = "#66A61E",
    "Risk type" = "#E69F00"
  )) +
  labs(
    title = "",
    subtitle = "",
    x = paste0("Dimension 1 (", dim1_1, "%)"),
    y = paste0("Dimension 2 (", dim2_1, "%)"),
    colour = "Variables"
  ) +
  theme_bw(base_size = 12) +
  theme(
    axis.title = element_text(face = "bold"),
    legend.position = "right",
    panel.grid.minor = element_blank()
  )

p_mca_1

ggsave(
  "MCA_categories_Prueba1_standardised_no_rare_no_reported_Benzecri_corrected.tiff",
  p_mca_1,
  width = 18,
  height = 16,
  units = "cm",
  dpi = figure_dpi,
  compression = figure_compression
)


# ============================================================
# 5. MCA ANALYSIS - MATRIX 2
# ============================================================

# Matrix 2 variables:
#   Specific mechanism: specific climate-related mechanism or stressor
#   Location: country or geographic region
#   Ecosystem: ecosystem type
#   Risk_Type: type of risk

data_2 <- read_excel(file_matrix_2, sheet = sheet_name)

names(data_2) <- clean_column_names(names(data_2))

data_mca_2 <- data_2 %>%
  select(
    `Specific mechanism`,
    Location,
    Ecosystem,
    Risk_Type
  ) %>%
  mutate(across(everything(), clean_text)) %>%
  filter(if_all(everything(), ~ !is.na(.) & . != ""))


# ------------------------------------------------------------
# 5.1 Standardise categories - Matrix 2
# ------------------------------------------------------------

data_mca_2 <- data_mca_2 %>%
  mutate(

    `Specific mechanism` = case_when(
      `Specific mechanism` %in% c(
        "Marine Heatwaves",
        "Marine heatwaves",
        "Marine heatwave"
      ) ~ "Marine heatwaves",

      `Specific mechanism` %in% c(
        "Sea level rise",
        "Sea-level rise",
        "SLR"
      ) ~ "Sea-level rise",

      `Specific mechanism` %in% c(
        "Ocean warming",
        "Ocean temperature increase",
        "Warming"
      ) ~ "Ocean warming",

      `Specific mechanism` %in% c(
        "Extreme precipitation",
        "Extreme rainfall",
        "Heavy rainfall"
      ) ~ "Extreme precipitation",

      `Specific mechanism` %in% c(
        "River discharge change",
        "River-discharge change",
        "River discharge variability"
      ) ~ "River discharge change",

      `Specific mechanism` %in% c(
        "Precipitation change",
        "Precipitation anomaly",
        "Precipitation decrease"
      ) ~ "Precipitation change",

      TRUE ~ `Specific mechanism`
    ),

    Location = case_when(
      Location %in% c("Brazil", "Brasil") ~ "Brazil",

      Location %in% c(
        "Colombia Pacific",
        "Colombia (Pacific)"
      ) ~ "Colombia (Pacific)",

      Location %in% c(
        "Colombia Caribbean",
        "Colombia (Caribbean)"
      ) ~ "Colombia (Caribbean)",

      Location %in% c(
        "Argentina Uruguay",
        "Argentina, Uruguay"
      ) ~ "Argentina, Uruguay",

      Location %in% c(
        "Argentina Brazil Uruguay",
        "Argentina, Brazil, Uruguay"
      ) ~ "Argentina, Brazil, Uruguay",

      TRUE ~ Location
    ),

    Ecosystem = case_when(
      Ecosystem %in% c("Coral reef", "Coral reefs") ~ "Coral reefs",
      Ecosystem %in% c("Mangrove", "Mangroves") ~ "Mangroves",
      Ecosystem %in% c("Kelp forest", "Kelp forests") ~ "Kelp forests",
      Ecosystem %in% c(
        "Shelf and pelagic",
        "Shelf/pelagic",
        "Shelf & pelagic"
      ) ~ "Shelf and pelagic",
      Ecosystem %in% c("Sandy beach", "Sandy beaches") ~ "Sandy beaches",
      Ecosystem %in% c("Coastal wetland", "Coastal wetlands") ~ "Coastal wetlands",
      Ecosystem %in% c("Tidal marsh", "Tidal marshes") ~ "Tidal marshes",
      TRUE ~ Ecosystem
    ),

    Risk_Type = case_when(
      Risk_Type %in% c(
        "Socio ecological",
        "Socio-ecological",
        "Socioecological"
      ) ~ "Socio-ecological",

      Risk_Type %in% c(
        "Biogeophysical",
        "Bio-geophysical"
      ) ~ "Biogeophysical",

      TRUE ~ Risk_Type
    )
  )


# ------------------------------------------------------------
# 5.2 Remove not reported values and rare categories - Matrix 2
# ------------------------------------------------------------

data_mca_2 <- data_mca_2 %>%
  filter(
    `Specific mechanism` != "Not reported",
    Location != "Not reported",
    Ecosystem != "Not reported",
    Risk_Type != "Not reported"
  )

if(remove_rare_categories){
  data_mca_2 <- data_mca_2 %>%
    mutate(across(
      everything(),
      ~ remove_rare_levels(., min_n = min_n)
    )) %>%
    drop_na()
}

data_mca_2 <- data_mca_2 %>%
  mutate(across(everything(), as.factor))

print(lapply(data_mca_2, table))
str(data_mca_2)
summary(data_mca_2)


# ------------------------------------------------------------
# 5.3 Run MCA - Matrix 2
# ------------------------------------------------------------

mca_2 <- MCA(data_mca_2, graph = FALSE)

corrected_inertia_2 <- mca_corrected_inertia(mca_2, data_mca_2)

print(corrected_inertia_2)


# ------------------------------------------------------------
# 5.4 Export results - Matrix 2
# ------------------------------------------------------------

export_mca_results(
  output_file = "MCA_results_Prueba2_standardised_no_rare_no_reported_corrected.xlsx",
  data_mca = data_mca_2,
  mca_res = mca_2,
  corrected_inertia = corrected_inertia_2
)


# ------------------------------------------------------------
# 5.5 Prepare plot data - Matrix 2
# ------------------------------------------------------------

coord_plot_2 <- as.data.frame(mca_2$var$coord)
names(coord_plot_2)[1:2] <- c("Dim1", "Dim2")
coord_plot_2$Category <- rownames(coord_plot_2)

map_variables_2 <- bind_rows(
  lapply(names(data_mca_2), function(v){
    data.frame(
      Category = levels(data_mca_2[[v]]),
      Variable = v
    )
  })
)

coord_plot_2 <- coord_plot_2 %>%
  left_join(map_variables_2, by = "Category") %>%
  mutate(
    Variable = recode(
      Variable,
      `Specific mechanism` = "Specific Climate Stressor",
      Location = "Location",
      Ecosystem = "Ecosystem",
      Risk_Type = "Risk type"
    )
  )

coord_plot_2$Variable <- factor(
  coord_plot_2$Variable,
  levels = c(
    "Specific Climate Stressor",
    "Location",
    "Ecosystem",
    "Risk type"
  )
)

table(coord_plot_2$Variable, useNA = "ifany")

dim1_2 <- round(corrected_inertia_2$Benzecri_percent[1], 1)
dim2_2 <- round(corrected_inertia_2$Benzecri_percent[2], 1)


# ------------------------------------------------------------
# 5.6 Plot and export figure - Matrix 2
# ------------------------------------------------------------

p_mca_2 <- ggplot(
  coord_plot_2,
  aes(x = Dim1, y = Dim2, colour = Variable)
) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.4) +
  geom_point(size = 1.8, alpha = 0.8) +
  geom_text_repel(
    aes(label = Category),
    size = 3,
    max.overlaps = Inf,
    box.padding = 0.25,
    point.padding = 0.15,
    segment.size = 0.2,
    show.legend = FALSE
  ) +
  scale_colour_manual(values = c(
    "Specific Climate Stressor" = "#6A5ACD",
    "Location" = "#4D4D4D",
    "Ecosystem" = "#66A61E",
    "Risk type" = "#E69F00"
  )) +
  labs(
    title = "",
    subtitle = "",
    x = paste0("Dimension 1 (", dim1_2, "%)"),
    y = paste0("Dimension 2 (", dim2_2, "%)"),
    colour = "Variables"
  ) +
  theme_bw(base_size = 12) +
  theme(
    axis.title = element_text(face = "bold"),
    legend.position = "right",
    panel.grid.minor = element_blank()
  )

p_mca_2

ggsave(
  "MCA_categories_Prueba2_standardised_no_rare_no_reported_Benzecri_corrected.tiff",
  p_mca_2,
  width = 18,
  height = 20,
  units = "cm",
  dpi = figure_dpi,
  compression = figure_compression
)


# ============================================================
# End of script
# ============================================================
